


library(mosaic) 
library(tidyverse)
library(jtools)
library(interplot)
library(stargazer)
require(gridExtra)
library(coefplot)
library(MASS)
library(modelsummary)

##### Import database and create the independent variable based on respondents' experimental group 

May23 <- Int_obl %>% mutate(treat1 = ifelse(FL_93_DO == "block_treatment" |   
                                            FL_95_DO == "block_treatment" | 
                                            FL_96_DO == "block_treatment" | 
                                            FL_97_DO == "block_treatment", 1, NA))

May23 <- May23 %>% mutate(treat2 = ifelse(FL_93_DO == "block_control" |   
                                            FL_95_DO == "block_control" | 
                                            FL_96_DO == "block_control" | 
                                            FL_97_DO == "block_control", 0, NA))

# sum the rows 

May23$treatment = rowSums(May23[,c(42, 43)], na.rm = TRUE)

# exclude those who haven't completed the survey 

May23 <- May23 %>% filter(dep1 != "NA")
May23 <- May23 %>% filter(dep3_tax != "NA")
May23 <- May23 %>% filter(dep2_spending != "NA")



##### create a composite dependent variable ##########

May23$dep_composite = rowSums(May23[,c(27, 28, 29)], na.rm = TRUE) 

# rescale the dep variables from 0 to 1

May23 <- May23 %>% mutate(dep_res = (dep_composite - min(dep_composite)) / (max(dep_composite) - min(dep_composite))) #dep_res = composite dependent variable
May23 <- May23 %>% mutate(dep1_res = (dep1 - min(dep1)) / (max(dep1) - min(dep1))) # dep1_res = no trade-off 
May23 <- May23 %>% mutate(dep2_res = (dep2_spending - min(dep2_spending)) / (max(dep2_spending) - min(dep2_spending))) # dep2_res = spending trade-off 
May23 <- May23 %>% mutate(dep3_res = (dep3_tax - min(dep3_tax)) / (max(dep3_tax) - min(dep3_tax))) # dep3_res = taxes trade-of

###################

# income from 1 to 10, from lower to higher income categories

May23$income[May23$income == 27] <- NA
May23$income[May23$income == 26] <- 1
May23$income[May23$income == 25] <- 2
May23$income[May23$income == 24] <- 3
May23$income[May23$income == 23] <- 4
May23$income[May23$income == 22] <- 5
May23$income[May23$income == 30] <- 6
May23$income[May23$income == 31] <- 7
May23$income[May23$income == 29] <- 8
May23$income[May23$income == 28] <- 9
May23$income[May23$income == 21] <- 10



####### HISTOGRAM ### Figure 3 main text ######

#distribution of the variable with ggplot


hist1 <- ggplot(May23, aes(x=dep1)) +
  geom_histogram(aes(y=..count../sum(..count..)), binwidth = 0.5) +
  labs(title ="Support for Debt Reduction (No Trade-Off)", 
       x = "", y = "Proportion of Observations") +
  scale_color_brewer(palette="Dark2") + 
  scale_x_continuous(breaks = c(0,1, 2, 3,4),label = c( "Totally disagree", "Disagree", "Neutral", "Agree", "Totally Agree")) +
  theme(plot.title = element_text(hjust = 0.5, size = 14))
hist1

hist2 <- ggplot(May23, aes(x=dep2_spending)) +
  geom_histogram(aes(y=..count../sum(..count..)), binwidth = 0.5) +
  labs(title ="Support for Debt Reduction (Spending Trade-Off)", 
       x = "The Government Should Reduce Public Debt", y = "") +
  scale_color_brewer(palette="Dark2") + 
  scale_x_continuous(breaks = c(0,1, 2, 3,4),label = c( "Totally disagree", "Disagree", "Neutral", "Agree", "Totally Agree")) +
  theme(plot.title = element_text(hjust = 0.5, size = 14))
hist2

hist3 <- ggplot(May23, aes(x=dep3_tax)) +
  geom_histogram(aes(y=..count../sum(..count..)), binwidth = 0.5) +
  labs(title ="Support for Debt Reduction (Taxes Trade-Off)", 
       x = "", y = "") +
  scale_color_brewer(palette="Dark2") + 
  scale_x_continuous(breaks = c(0,1, 2, 3,4),label = c( "Totally disagree", "Disagree", "Neutral", "Agree", "Totally Agree")) +
  theme(plot.title = element_text(hjust = 0.5, size = 14))
hist3


grid.arrange(hist1, hist2, hist3, ncol=3)

###########################################################################


####### testing H1

## with quotas ##### TABLE B1 Appendix

m0 <-lm(dep_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(m0) 


m1 <-lm(dep1_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(m1)


m2 <-lm(dep2_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(m2)

m3 <-lm(dep3_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(m3)



########### COEFFICIENT PLOT ##################### FIGURE 4 main text #######################

multiplot(m0,m1,m2,m3, coefficients="treatment", title = "",   sort = "alphabetical", 
          xlab = "", ylab = "", plot.shapes = TRUE, 
          innerCI = 0,
          names = c("Composite variable", "No trade-off","Spending trade-off", "Taxes trade-off")) + 
  theme_bw() +
  theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(), legend.title = element_text(size=14),
        legend.text = element_text(size=13)) +
  scale_y_discrete(name="", labels=" International treatment") + 
  guides(shape = guide_legend(reverse = TRUE), color = "none") +
  xlim(-0.08, 0.09) + 
  theme(axis.text.y = element_text(size=14)) +
 scale_colour_manual(values = c("black", "black", "black", "black"))



#################################################################################


#################### testing H2

#recode the variable European_countries so that higher values of the moderator correspond to higher trust

May23$European_countries[May23$European_countries == 1] <- 3 
May23$European_countries[May23$European_countries == 2] <- 2 
May23$European_countries[May23$European_countries == 4] <- 1 
May23$European_countries[May23$European_countries == 5] <- 0 


May23 <- May23 %>% mutate(EU_res = (European_countries - min(European_countries)) / (max(European_countries) - min(European_countries)))

######### TABLE B2 Appendix ########## with quotas

EU0 <-lm(dep_res ~ EU_res + age + as.factor(sex) + as.factor(residence) + as.factor(education),
                data = May23)
summ(EU0)

EU1 <-lm(dep1_res ~ EU_res + age + as.factor(sex) + as.factor(residence) + as.factor(education), 
                data = May23)
summ(EU1)

EU2 <-lm(dep2_res ~ EU_res + age + as.factor(sex) + as.factor(residence) + as.factor(education),
                data = May23)
summ(EU2)

EU3 <-lm(dep3_res ~ EU_res + age + as.factor(sex) + as.factor(residence) + as.factor(education),
                data = May23)
summ(EU3)


## INTERACTION regression ##### TABLE B3 Appendix ########

m_EUint <-lm(dep_res ~ treatment*EU_res + age + as.factor(residence)
             + as.factor(education) + as.factor(sex) , data = May23)
summary(m_EUint)
summ(m_EUint)


m1_EUint <-lm(dep1_res ~ treatment*EU_res + age + as.factor(residence) 
              + as.factor(education) + as.factor(sex) , data = May23)
summary(m1_EUint)
summ(m1_EUint)


m2_EUint<-lm(dep2_res ~ treatment*EU_res + age + as.factor(residence) 
             + as.factor(education) + as.factor(sex) , data = May23)
summary(m2_EUint)
summ(m2_EUint)


m3_EUint <-lm(dep3_res ~ treatment*EU_res + age + as.factor(residence) 
              + as.factor(education) + as.factor(sex) , data = May23)
summary(m3_EUint)
summ(m3_EUint)



######## Figure 5 main text ##############################

interplot(m_EUint, "treatment", "EU_res", ci = 0.90 ) +
  xlab("View of the Governments in Other European Countries") +
  ylab("Estimated Coefficient of Treatment") +
  # Change the background
  theme_bw() +
  # Add the title
  
  theme(plot.title = element_text(face="bold")) +
  
  theme(plot.title = element_text(hjust = 0.5, size = 16), 
        axis.title.x = element_text(hjust = 0.5, size = 14), 
        axis.title.y = element_text(hjust = 0.5, size = 14)) +
  scale_x_continuous(breaks = c(0,0.333, 0.666, 1), label = c( "Very Unfavorable ", "Somewhat Unfavorable", "Somewhat Favorable", "Very Favorable")
)


###########################################################################################�


#mechanisms ########## Table B4 Appendix #############


# (1) 
mech_4<-lm(moral_duty_1 ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(mech_4) # logic of appropriateness 

# (2)
mech_2<-lm(reciprocity_1 ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(mech_2) #reciprocity


#(3)
mech_5<-lm(signalling_1 ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(mech_5) # signalling

# (4) 
mech_1<-lm(efficacy_1 ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(mech_1) # effectiveness

#(5)
mech_3<-lm(credibility_1 ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(mech_3) #logic of consequences





### coef plot ######### FIGURE 6 main text #######################

multiplot(mech_4,mech_2,mech_3, mech_1 ,mech_5,  coefficients="treatment", innerCI = 0, title = "",   sort = "alphabetical",
          xlab = "", ylab = "", plot.shapes = TRUE, 
          names = c("Moral duty", "Reciprocity","Credibility", "Effectiveness", "Signaling")) + theme_bw() +
  theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(), legend.title = element_text(size=14),
        legend.text = element_text(size=13)) +
  scale_y_discrete(name="", labels=" International treatment") + 
  theme(axis.text.y = element_text(size=14)) +
  scale_colour_manual(values = c("black", "black", "black", "black", "black"))



#####

mean(May23$moral_duty_1, na.rm = TRUE)
mean(May23$efficacy_1, na.rm = TRUE)
mean(May23$signalling_1, na.rm = TRUE)
mean(May23$credibility_1, na.rm = TRUE)
mean(May23$reciprocity_1, na.rm = TRUE)


################## PRIOR INFO ##############  (pre-registered)


#recode the variable so that higher values of the moderator correspond to higher knowledge (originally from 4 to 1)

May23$prior_info[May23$prior_info == 4] <- 0 
May23$prior_info[May23$prior_info == 3] <- 20 # 
May23$prior_info[May23$prior_info == 2] <- 2 
May23$prior_info[May23$prior_info == 1] <- 3 
May23$prior_info[May23$prior_info == 20] <- 1 

May23 <- May23 %>% mutate(priorinfo_res = (prior_info - min(prior_info)) / (max(prior_info) - min(prior_info)))


#preliminary

info_pre <-lm(dep_res ~ priorinfo_res + age + as.factor(sex) + as.factor(residence) + as.factor(education) , data = May23)
summ(info_pre)

#interaction ############ TABLE D1 Appendix ############

m0_INFOint <-lm(dep_res ~ priorinfo_res*treatment , data = May23)
summary(m0_INFOint)
summ(m0_INFOint)


m1_INFOint <-lm(dep1_res ~ treatment*priorinfo_res, data = May23)
summary(m1_INFOint)
summ(m1_INFOint)


m2_INFOint<-lm(dep2_res ~ treatment*priorinfo_res, data = May23)
summary(m2_INFOint)
summ(m2_INFOint)


m3_INFOint <-lm(dep3_res ~ treatment*priorinfo_res, data = May23)
summary(m3_INFOint)
summ(m3_INFOint)




######## alternative explanation ########### DISCUSSION SECTION ###################�

# support rule : originally from 0 to 4 support decreasing  
#flip the variable (variable increases, support increases)

May23$support_rule[May23$support_rule == 0] <- 20
May23$support_rule[May23$support_rule == 1] <- 21
May23$support_rule[May23$support_rule == 2] <- 22
May23$support_rule[May23$support_rule == 3] <- 23
May23$support_rule[May23$support_rule == 4] <- 24

May23$support_rule[May23$support_rule == 20] <- 4
May23$support_rule[May23$support_rule == 21] <- 3
May23$support_rule[May23$support_rule == 22] <- 2
May23$support_rule[May23$support_rule == 23] <- 1
May23$support_rule[May23$support_rule == 24] <- 0


May23_s <- May23 %>% filter(support_rule != "NA")
May23_s <- May23_s %>% mutate(supp_res = (support_rule - min(support_rule)) / (max(support_rule) - min(support_rule)))

### Table D2 Appendix 

m_s0 <-lm(supp_res ~ treatment, data = May23_s)
summ(m_s0)

m_s1 <-lm(supp_res ~ treatment + EU_res + age + as.factor(sex) + income + as.factor(education) + as.factor(residence) 
          + right_left, data = May23_s)
summ(m_s1)

###### FIGURE D1 Appendix ######

hist4 <- ggplot(May23_s, aes(x=support_rule)) +
  geom_histogram(aes(y=..count../sum(..count..)), binwidth = 0.5) +
  labs(title ="Support for the New Fiscal Rule", 
       x = "", y = "") +
  scale_color_brewer(palette="Dark2") + 
  scale_x_continuous(breaks = c(0,1, 2, 3,4),label = c( "Totally against", "Against", "Neutral", "In favour", "Totally in favour")) +
  theme(plot.title = element_text(hjust = 0.5, size = 14))
hist4

######## Table D3 Appendix #########

# legitimacy 0 to 4, increasing support for European countries deciding for Italy 

May23_l <- May23_s %>% filter(legitimacy != "NA")

May23_l <- May23_l %>% mutate(leg_res = (legitimacy - min(legitimacy)) / (max(legitimacy) - min(legitimacy)))


m_l0 <-lm(leg_res ~ treatment, data = May23_l)
summ(m_l0)

m_l1 <-lm(leg_res ~ treatment + EU_res + age + as.factor(sex) + income + as.factor(education) + as.factor(residence) 
          + right_left, data = May23_l)
summ(m_l1)
summary(m_l1)


################# BALANCE TABLE, table A1 Appendix ########## 

library(tables)
library(modelsummary)


May23_balance <- May23 %>%
  mutate(Age = age,
         Gender = sex,
         Education = education,
         "Region of residence" = residence,
         "Left-Right orientation" = right_left,
         "View of other European countries" = EU_res,
         "Prior fiscal knowledge" = priorinfo_res,
         Treatment = ifelse(treatment == 1, "Treatment", "Control"))



May23_balance  <- May23_balance %>% subset(select = c(52:59))


balance_table <- datasummary_balance(~ Treatment,
                                     data = May23_balance)

balance_table


############## EQUIVALENCE TESTING #### Tables from E1 to E4 Appendix 

library(TOSTER)

May_a <- May23 %>% filter(treatment == 0)
May_b <- May23 %>% filter(treatment == 1)


# composite var

m1 <- mean(May_a$dep_res, na.rm = TRUE) 
m2 <- mean(May_b$dep_res, na.rm = TRUE) 


s1 <- sd(May_a$dep_res, na.rm = TRUE) 
s2 <- sd(May_b$dep_res, na.rm = TRUE) 

TOSTER::tsum_TOST(m1 = m1, 
                  m2 = m2, 
                  sd1 = s1, 
                  sd2 = s2,
                  n1 = 616, 
                  n2 = 602, 
                  low_eqbound = -0.05, 
                  high_eqbound = 0.05)



# No trade-off

m1 <- mean(May_a$dep1_res, na.rm = TRUE) 
m2 <- mean(May_b$dep1_res, na.rm = TRUE) 


s1 <- sd(May_a$dep1_res, na.rm = TRUE) 
s2 <- sd(May_b$dep1_res, na.rm = TRUE) 

TOSTER::tsum_TOST(m1 = m1, 
                  m2 = m2, 
                  sd1 = s1, 
                  sd2 = s2,
                  n1 = 616, 
                  n2 = 602, 
                  low_eqbound = -0.05, 
                  high_eqbound = 0.05)

#### spending

m1 <- mean(May_a$dep2_res, na.rm = TRUE)
m2 <- mean(May_b$dep2_res, na.rm = TRUE)



s1 <- sd(May_a$dep2_res, na.rm = TRUE) 
s2  <- sd(May_b$dep2_res, na.rm = TRUE) 




TOSTER::tsum_TOST(m1 = m1, 
                  m2 = m2, 
                  sd1 = s1, 
                  sd2 = s2,
                  n1 = 616, 
                  n2 = 602, 
                  low_eqbound = -0.05, 
                  high_eqbound = 0.05)

# taxes

m1 <- mean(May_a$dep3_res, na.rm = TRUE) #1.313312
m2 <- mean(May_b$dep3_res, na.rm = TRUE) # 1.393688


s1<- sd(May_a$dep3_res, na.rm = TRUE) # 1.071835
s2  <- sd(May_b$dep3_res, na.rm = TRUE) # 1.127257



TOSTER::tsum_TOST(m1 = m1, 
                  m2 = m2, 
                  sd1 = s1, 
                  sd2 = s2,
                  n1 = 616, 
                  n2 = 602, 
                  low_eqbound = -0.05, 
                  high_eqbound = 0.05)


##### Attention check ############ Appendix F ############�


## 1) Table F1, duration (exclude outliers)

colnames(May23)[1] <- "duration"

mean_duration <- mean(May23$duration)
sd_duration <- sd(May23$duration)

May23 <- May23 %>% mutate(standardized_duration = ((duration - mean_duration)/ sd_duration))

# mean is 0 and standard deviation is 1,approximately 95% of the data falls within 1.96 standard deviations from the mean. 
# Values beyond this threshold in either direction can be considered as potentially unusual or extreme.

May23_att <- May23 %>% filter(standardized_duration >= -1.96 & standardized_duration <= 1.96)

m0_att <-lm(dep_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att)
summ(m0_att) 


m1_att <-lm(dep1_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att)
summ(m1_att)


m2_att <-lm(dep2_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att)
summ(m2_att)


m3_att <-lm(dep3_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att)
summ(m3_att)


### 2) duration vignette ########### Table F2 ## 


May23_att_vignette <- May23 %>% rowwise() %>%  mutate(duration_vignette = sum(`delay_next2_Page Submit`, `delay_next_t2_Page Submit`, na.rm= TRUE))
May23_att_vignette <- May23_att_vignette %>% filter(duration_vignette > 14) # exclude 10% quickest 


m0_v <-lm(dep_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att_vignette)
summ(m0_v) 


m1_v <-lm(dep1_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att_vignette)
summ(m1_v)


m2_v <-lm(dep2_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att_vignette)
summ(m2_v)


m3_v <-lm(dep3_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att_vignette)
summ(m3_v)


## 3)  varible name: att_check_1 (answer "2" to the next question)

May23_att3 <- May23 %>% filter(att_check_1 == 2) #971 respondents didn't fail the attention check 

####### TABLE F3 Appendix with quotas

m0_att3 <-lm(dep_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att3)
summ(m0_att3) 


m1_att3 <-lm(dep1_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att3)
summ(m1_att3)


m2_att3 <-lm(dep2_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att3)
summ(m2_att3)


m3_att <-lm(dep3_res ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23_att3)
summ(m3_att)
summary(m3_att)

sd(May23_att3$dep3_res)


#################### Other moderators ##########################################
########## APPENDIX G #########################

### right_left: from 0 (left) to 10 (right) Table G1


m_lr <-lm(dep_res ~ treatment*right_left + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(m_lr)


m1_lr <-lm(dep1_res ~ treatment*right_left + age + as.factor(sex) + as.factor(residence) + as.factor(education) , data = May23)
summ(m1_lr)


m2_lr <-lm(dep2_res ~ treatment*right_left + age + as.factor(sex) + as.factor(residence) + as.factor(education) , data = May23)
summ(m2_lr)


m3_lr <-lm(dep3_res ~ treatment*right_left + age + as.factor(sex) + as.factor(residence) + as.factor(education) , data = May23)
summ(m3_lr)


##### education Table G2

m_edu <-lm(dep_res ~ treatment*as.factor(education) + age + as.factor(sex) + as.factor(residence), data = May23)
summ(m_edu)


m1_edu <-lm(dep1_res ~ treatment*as.factor(education) + age + as.factor(sex) + as.factor(residence) , data = May23)
summ(m1_edu)


m2_edu <-lm(dep2_res ~ treatment*as.factor(education) + age + as.factor(sex) + as.factor(residence)  , data = May23)
summ(m2_edu)


m3_edu <-lm(dep3_res ~ treatment*as.factor(education) + age + as.factor(sex) + as.factor(residence)  , data = May23)
summ(m3_edu)

#################### income ############ Table G3
 

May23$income <- as.numeric(May23$income)

m_inc <-lm(dep_res ~ treatment*income + age + as.factor(sex) + as.factor(residence) + as.factor(education), data = May23)
summ(m_inc)


m1_inc <-lm(dep1_res ~ treatment*income + age + as.factor(sex) + as.factor(residence) + as.factor(education) , data = May23)
summ(m1_inc)


m2_inc <-lm(dep2_res ~ treatment*income + age + as.factor(sex) + as.factor(residence) + as.factor(education) , data = May23)
summ(m2_inc)


m3_inc <-lm(dep3_res ~ treatment*income + age + as.factor(sex) + as.factor(residence)+ as.factor(education)  , data = May23)
summ(m3_inc)



####### Ordinal Logistic Regression ##### Appendix H #####

#Table H1

m0_logistic <-polr(as.factor(dep_res) ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), 
                   data = May23, Hess=TRUE)
modelsummary(m0_logistic, stars=TRUE)


m1_logistic <-polr(as.factor(dep1_res) ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), 
                   data = May23, Hess=TRUE)
modelsummary(m1_logistic, stars=TRUE)


m2_logistic <-polr(as.factor(dep2_res) ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), 
                   data = May23, Hess=TRUE)
modelsummary(m2_logistic, stars=TRUE)

m3_logistic <-polr(as.factor(dep3_res) ~ treatment + age + as.factor(sex) + as.factor(residence) + as.factor(education), 
                   data = May23, Hess=TRUE)
modelsummary(m3_logistic, stars=TRUE)


# Table H2 (interaction)

m_EUint_logistic <-polr(as.factor(dep_res) ~ treatment*EU_res +  age + as.factor(residence) + as.factor(education) + 
                          as.factor(sex),
                         data = May23, Hess=TRUE)

modelsummary(m_EUint_logistic, stars=TRUE)


m1_EUint_logistic <-polr(as.factor(dep1_res) ~ treatment*EU_res +  age + as.factor(residence) + as.factor(education) + 
           as.factor(sex) , data = May23, Hess=TRUE)
modelsummary(m1_EUint_logistic, stars=TRUE)



m2_EUint_logistic <-polr(as.factor(dep2_res) ~ treatment*EU_res +  age + as.factor(residence) + as.factor(education) + 
                           as.factor(sex), data = May23, Hess=TRUE)
modelsummary(m2_EUint_logistic, stars=TRUE)



m3_EUint_logistic <-polr(as.factor(dep3_res) ~ treatment*EU_res +  age + as.factor(residence) + as.factor(education) + 
                           as.factor(sex), data = May23, Hess=TRUE)
modelsummary(m3_EUint_logistic, stars=TRUE)


# Table H3 prior information

m0_INFOint <-polr(as.factor(dep_res) ~ priorinfo_res*treatment , data = May23)
modelsummary(m0_INFOint, stars=TRUE)

m1_INFOint <-polr(as.factor(dep1_res) ~ treatment*priorinfo_res, data = May23)
modelsummary(m1_INFOint, stars=TRUE)


m2_INFOint<-polr(as.factor(dep2_res) ~ treatment*priorinfo_res, data = May23)
modelsummary(m2_INFOint, stars=TRUE)

m3_INFOint <-polr(as.factor(dep3_res) ~ treatment*priorinfo_res, data = May23)
modelsummary(m3_INFOint, stars=TRUE)

# Table H4 support for the rule 

m_s0_logistic  <-polr(as.factor(supp_res) ~ treatment, data = May23_s, Hess=TRUE)
modelsummary(m_s0_logistic, stars=TRUE)

m_s1_logistic  <-polr(as.factor(supp_res) ~ treatment + EU_res + age + as.factor(sex) + income + as.factor(education) + as.factor(residence) 
          + right_left, data = May23_s, Hess=TRUE)
modelsummary(m_s1_logistic, stars=TRUE)


# Table H5 perceived legitimacy

m_l0_logistic <-polr(as.factor(leg_res) ~ treatment, data = May23_l, Hess=TRUE)
modelsummary(m_l0_logistic, stars=TRUE)

m_l1_logistic <-polr(as.factor(leg_res) ~ treatment + EU_res + age + as.factor(sex) + income + as.factor(education) + as.factor(residence) 
          + right_left, data = May23_l, Hess=TRUE)
modelsummary(m_l1_logistic, stars=TRUE)



## census Appendix C Table C1 

#female

num_females <- May23 %>% 
  filter(sex == 2) %>%
  summarise(count = n())

print(num_females) # 616

percent_female <- 616*100/1218 #n. female*100 / sample size

# age

May23$age <- as.numeric(May23$age)

age_groups <- May23 %>%
  mutate(age_group = case_when(
    age >= 18 & age <= 29 ~ "18-29",
    age >= 30 & age <= 39 ~ "30-39",
    age >= 40 & age <= 49 ~ "40-49",
    age >= 50 & age <= 59 ~ "50-59",
    age >= 60 ~ "60+",
    TRUE ~ NA_character_)) %>%
  group_by(age_group) %>%
  summarise(
    count = n(),
    percentage = (n() / nrow(May23)) * 100
  ) %>%
  filter(!is.na(age_group))  # Exclude any NA groups

print(age_groups)

# region 

May23$residence

#1 Abruzzo
#2 Basilicata
#3 Calabria
#4 Campania
#5 Emilia-Romagna
#6 Friuli Venezia Giulia
#7 Lazio
#8 Liguria
#9 Lombardia
#10 Marche
#11 Molise
#12 Piemonte
#14 Puglia
#15 Sardegna
#16 Sicilia
#17 Toscana
#18 Trentino-Alto Adige
#19 Umbria
#20 Valle d'Aosta
#21 Veneto

May23 <- May23 %>%
  mutate(macroregion = case_when(
    residence %in% c(1, 2, 3, 4, 10, 14) ~ "Sud",
    residence %in% c(5, 6, 18, 21) ~ "Nord-Est",
    residence %in% c(7, 10, 17, 19) ~ "Centro",
    residence %in% c(15, 16) ~ "Isole",
    residence %in% c(8, 9, 12, 20) ~ "Nord-Ovest",
    TRUE ~ "Unknown"
  ))

macroregion_summary <- May23 %>%
  group_by(macroregion) %>%
  summarise(
    count = n(),
    percentage = (n() / nrow(May23)) * 100
  ) %>%
  filter(macroregion != "Unknown")  # Exclude any unknown macroregions

print(macroregion_summary)
